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A quintessence scalar field or cosmon interacting with neutrinos can have important effects on 
cosmological structure formation. Within growing neutrino models the coupling becomes effective 
only in recent times, when neutrinos become nonrelativistic, stopping the evolution of the cosmon. 
This can explain why dark energy dominates the Universe only in a rather recent epoch by relating 
the present dark energy density to the small mass of neutrinos. Such models predict the presence 
of stable neutrino lumps at supercluster scales (~ 200 Mpc and bigger), caused by an attractive 
force between neutrinos which is stronger than gravity and mediated by the cosmon. We present a 
method to follow the initial nonlinear formation of neutrino lumps in physical space, by integrating 
numerically on a 3D grid nonlinear evolution equations, until virialization naturally occurs. As a 
first application, we show results for cosmologies with final large neutrino average mass ~ 2 eV: 
in this case, neutrino lumps indeed form and mimic very large cold dark matter structures, with 
a typical gravitational potential 10~ 5 for a lump size ~ 10 Mpc, and reaching larger values for 
lumps of about 200 Mpc. A rough estimate of the cosmological gravitational potential at small k 
in the nonlinear regime, = 10~ 6 (k/k )~ 2 , 1.2 ■ 10~ 2 h/Mpc < ko < 7.8 ■ 10 -2 h/Mpc, turns 
out to be many orders of magnitude smaller than an extrapolation of the linear evolution of density 
fluctuations. The size of the neutrino-induced gravitational potential could modify the spectrum of 
CMB anisotropics for small angular momenta. 



I. INTRODUCTION 

The presence of an interaction between a quintessence 
scalar field or cosmon and other species in the Universe 
[H Q influences the nature and properties of dark en- 
ergy, with relevant effects on structure formation 
Recently, iV-body simulations have been performed for 
dark matter particles interacting with dark energy 

tWe concentrate here on growing neutrino models 
19], where the ncutrino-cosmon coupling can explain 
the "why now?" problem of dark energy. It has been re- 
cently shown [2(| that this predicts the existence of very 
large (supercluster) structures. 

Indeed, the key ingredient of growing neutrino 
quintessence is the presence of a coupling between dark 
energy and neutrinos. The latter have a mass that grows 
with the evolution of the Universe - m u {4>) being a func- 
tion of the cosmon field (p. The cosmon-neutrino coupling 
P is given by the logarithmic derivative j3 — —din m v /d(j). 
Since neutrino masses become cosmologically relevant 
only for redshift z«5, this framework can naturally an- 
swer the question why only in recent times dark energy 
leads the Universe expansion to accelerate, thus provid- 
ing a solution to the coincidence problem. 

In these models, as long as neutrinos stay relativis- 
tic, the coupling plays no role and dark energy tracks 
the background, along the attractor trajectories char- 
acterizing the cosmon evolution in the presence of an 
exponential potential 0, EH-SH. When neutrinos be- 
come nonrelativistic, the coupling between neutrinos and 
quintessence becomes relevant and almost stops the evo- 



lution of the cosmon. Then dark energy starts to resem- 
ble a cosmological constant with roughly the value of the 
exponential potential at the end of the attractor era. The 
transition from the attractor solution to an almost static 
solution is therefore strictly connected to a cosmological 
event, that is neutrinos becoming nonrelativistic. This 
"trigger event" leads dark energy to dominate over cold 
dark matter, naturally starting the recent era of acceler- 
ated expansion. 

Growing neutrino quintessence requires a cosmon- 
neutrino interaction somewhat stronger than gravity - 
typically the attraction between nonrelativistic neutri- 
nos exceeds gravity by a factor 10 3 . This coupling is 
substantially larger than a possible cosmon coupling to 
atoms. This may be motivated by the particular parti- 
cle physics mechanism responsible for the neutrino mass, 
which typically involves a heavy singlet field and not only 
the standard Higgs doublet [19(. Even much larger neu- 
trino couplings leading to a strongly coupled "acceleron- 
neutrino fluid" have been investigated within mass vary- 
ing neutrino models 0, [25 , 26 1 . In particular, mass vary- 
ing neuitrino models employ a scalar field with a mass 
much larger than the Hubble parameter. For growing 
neutrino models, in contrast, the time dependent cosmon 
mass equals the Hubble parameter up to a factor of or- 
der one, similar to many models of coupled quintessence. 
(For coupled quintessence with neutrinos at the linear 
level see ref. [Io|, 20|, 27-29[.) For this reason, the cosmon 
and the neutrinos always have to be treated as separate 
ingredients rather than as a common fluid. We remark 
that the small time dependent cosmon mass follows nat- 
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urally from possible explanations of an exponential po- 
tential in forms of an asymptotically vanishing dilatation 
anomaly (30j . It requires no additional fine tuning of 
particle physics parameters. 

Furthermore, in growing neutrino cosmologies the cou- 
pling is ineffective for most of the cosmological evolution 
and only becomes active when neutrinos become non- 
relativistic, relating naturally dark energy and neutrino 
properties. In view of bounds on the present neutrino 
mass, m„(to) < 2.3eV [3l|, and the time dependence of 
m u , which makes the mass even smaller in the past, the 
time when neutrinos become nonrelativistic is typically 
in the recent history of the Universe, say zj\jr ~ 5 — 10 
18]. It is only from this time on that neutrinos start feel- 
ing the effects of the coupling, manifesting effectively as a 
new attractive interaction between neutrinos. The 'fifth 
force' responsible for the formation of neutrino lumps 
"switches on" only in rather recent cosmology. 

Neutrino fluctuations on length scales larger than the 
free streaming length are still present at Znr,) and they 
start growing for z < znr with a large growth rate. As 
illustrated in j2fj, this opens the possibility that neu- 
trino perturbations rapidly grow nonlinear on superclus- 
ter scales and beyond. Neutrino nonlinear fluctuations 
later turn into bound neutrino lumps of the type dis- 
cussed in 32], thus opening a window for observable ef- 
fects of the growing neutrino scenario. The linear analy- 
sis [13] has already provided an estimate of this effect as a 
function of redshift and scale. Typically large scale neu- 
trino fluctuations with size £ 10 Mpc become nonlinear 
at a redshift z ~ 1. 

As noted in [2(J, a continuation of the linear evolu- 
tion beyond the time when the neutrino fluctuations are 
of order unity can easily produce erroneous results. In 
the linear approximation, the neutrino density contrast 
would quickly reach huge values, producing a very large 
gravitational potential. Then a very strong ISW effect 
would seemingly indicate a strong conflict with the cos- 
mic microwave anisotropies (cf. ref. [29| for a linear 
analysis). The true physics differs very strongly from the 
linear behavior, for example by the asymmetry between 
very large positive density contrasts, while a negative 
density contrast is bounded by 8 V > —1 since the neu- 
trino density cannot be negative. In order to understand 
the true gravitational potentials that will be generated 
by the neutrino lumps, one has to understand the non- 
linear dynamics of how local neutrino lumps form, how 
they are distributed in mass and size and how they may 
merge into larger lumps as time goes on. 

In this work we investigate the formation of individual 
neutrino lumps at a nonlinear level and in the Newtonian 
limit. One may wonder if a spherical collapse approach 
suffices to give a meaningful description of the lumps. We 
find that this is not the case, as the major force driving 
neutrinos to collapse is not gravity but the additional 



fifth force introduced by the coupling to the cosmon and 
dominant once neutrinos decouple from the background 
expansion. The evolution of the effective cosmic scale 
factor for the space occupied by the neutrino lump is only 
a subleading effect, in contrast to the formation of dark 
matter halos. We have therefore developed a numerical 
method for solving the hydrodynamic equations for the 
neutrino fluid, coupled to the cosmon and gravity. We 
have also included dark matter, but this is a subleading 
effect. 

The nonlinear analysis developed in this work provides 
a self-consistent way of analyzing the growth of neutrino 
perturbations in growing neutrino models. We show that 
neutrinos indeed form stable structures on large sub- 
horizon scales and we estimate the properties of the lump 
as a function of redshift and scale. In particular, we com- 
pute the gravitational potential of the lump at a time 
when the collapse ends due to virialization. This is a key 
quantity for an estimate of the effects of neutrino lumps 
on large scale structure - as large scale peculiar velocities 
- or on the cosmic microwave anisotropies in form of the 
integrated Sachs - Wolfe (ISW) effect. 

This paper is organized as follows. In section |TT] we 
recall the framework of growing neutrino cosmologies in 
which neutrino lumps form. In section IIIII we introduce 
the set of equations describing the evolution of neutrino 
ovcrdcnsities at a nonlinear level. We comment on the 
methods used in the numerical integration and present 
our results for the case of large present neutrino masses 
in section IIVI Section IVII discusses the initial conditions 
used for the nonlinear analysis. In section fVTTl we relate 
the nonlinear equations to relativistic linear ones. Finally 
we draw our conclusions in section IVIIII 

II. GROWING NEUTRINO COSMOLOGIES 

Growing neutrino models are described by the set of 
equations illustrated in 20j both for the evolution of the 
homogeneous and isotropic background and for linear 
perturbations. Here we recall for convenience the es- 
sential ingredients characterizing these models. At the 
background level, the Universe evolves in time according 
to the Friedmann and acceleration equations: 

*-(7)"-t5>-? w 

x / a 

and 

a 6 ' 

a 

where primes denote derivatives with respect to confor- 
mal time r, the sum is taken over all components a of the 
energy density in the Universe. We use k — for a spa- 
tially flat background. The equation of state w a is related 
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to the energy density p a for each species in the usual way, 
w a = p a l p a - A crucial ingredient in this model is the de- 
pendence of the neutrino mass on the cosmon field <f>, as 
encoded in the dimensionless cosmon-neutrino coupling 

P, 



d In to,. 



(3) 



For increasing 
with time 



and j3 < the neutrino mass increases 



(4) 



where rh v is a constant. The coupling /3 is chosen here to 
be a constant but can be, in general, a function of 0, as 
proposed in [l9[ within a particle physics model, leading 
to similar effects. The cosmon field <fi is normalized in 
units of the reduced Planck mass M = {SitGn)- 1 / 2 , and 
/3 2 gives the strength of the cosmon mediated interaction. 
The case f3 ~ 1 corresponds to a strength comparable to 
gravity. For a given cosmological model with a given time 
dependence of <f>, one can determine the time dependence 
of the neutrino mass to„ (t) . For three degenerate neutri- 
nos the present value of the neutrino mass m u (to) can be 
related to the energy fraction in neutrinos (/( « 0.72) 



94 eVh 2 



(5) 



The dynamics of the cosmon can be inferred from the 
Klein Gordon equation, now including an extra source 
due to the neutrino coupling, 



2H4>' + a' 



dU 



a 2 /3(Pv - 3p„) 



(G) 



with p u and p„ = w„p„ the energy density and pressure 
of the neutrinos. We choose an exponential potential 

13, [H-dl m : 



The sum of the energy momentum tensors for neutri- 
nos and dark energy is conserved, but not the separate 
parts. We neglect a possible cosmon coupling to Cold 
Dark Matter (CDM), so that p' c = -3Hp c . 

Given the potential ([7]), the evolution equations for the 
different species can be numerically integrated, providing 
the background evolution shown in figH] (Here we choose 
P = —52, a = 10 as in the original proposal [l8]). The 
initial pattern is a typical early dark energy model, since 
neutrinos are still relativistic and almost massless, with 
Vv = Pu/3 so that the coupling term in eq.® and (0 
[TO")) vanishes. Dark energy is still subdominant and falls 
into the attractor provided by the exponential potential 
(see [HSim for details), in which it tracks the dominant 
background component with an early dark energy frac- 
tion Qh = n/a 2 and n = 3(4) for the matter (radiation) 
dominated era. Radiation dominates until matter radia- 
tion equality, then CDM takes over. As the mass of the 
neutrinos increases with time, the coupling term ~ f3p v 
in the evolution equation for the cosmon ([6]) (or equiva- 
lently in ©[Til])) starts to play a significant role, kicking cf) 
out of the attractor as soon as neutrinos become nonrel- 
ativistic. In figUJthis is visible in the modified behavior 
of p v and pcf, for z < 10. Subsequently, small decaying 
oscillations characterize the <j) — v coupled fluid and the 
two components reach almost constant values. The va- 
lues of the energy densities today are in agreement with 
observations, once the precise crossing time for the end 
of the scaling solution has been fixed by an appropriate 
choice of the coupling /3. At present the neutrinos are 
still subdominant with respect to CDM, though in the 
future they will take the lead. 

For completeness, note that the unperturbed neutrino 
pressure reads 
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l 2 dqdQ 



fo(q) 



(11) 



V(<f>) = M 2 U{(j>) 



M 4 £ 



-OL(f> 



(7) 



where the constant a is one of the free parameters of our 
model and determines the amount of nonnegligible dark 
energy at early times. Current bounds constrain it to be 
of the order a ~ 10 or bigger [34 1. 

The homogeneous energy density and pressure of the 
scalar field 4> are defined in the usual way as 



qh is the comoving 3-momentum, 
,,(</>) 2 a 2 , / is the phase space distri- 



where q = ap = 

e = e O) = + n ' 

bution and fo its zeroth-order term (Fermi-Dirac distri- 
bution) . The neutrino energy density can either be given 
by solving the conservation equation (|10[) or equivalently 
via the integral: 



(12) 



p v = a q dqdQe(ct>)f a (q) 



P<f = ^o+V(<P) , P* = £z-VW , w^ = ^ . (8) 



2a 



2a 



Finally, we can express the conservation equations for 
dark energy and growing neutrinos as follows 0, Q : 



Pu 



-mil + w^pt+Wtl-Zw^pv , (9) 

-m{i + w v ) Pu - (3<j)\i-iw v )p u . (io) 



III. NONLINEAR EVOLUTION EQUATIONS 

Once neutrinos become nonrelativistic {z ~ 5 — 10), 
they start feeling an attractive force stronger than grav- 
ity and mediated by the cosmon field. Neutrino pertur- 
bations rapidly grow and become nonlinear at a redshift 
z ~ 1 — 2 [20l | . when they might form stable lumps, whose 
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FIG. 1. Energy densities of neutrinos (solid, red), cold dark 
matter (long-dashed, green), dark energy (dot-dashed, blue) 
and photons (short dashed, black) are plotted vs redshift. For 
all plots we take a constant ft — —52, with a — 10 and large 
average neutrino mass m v = 2.11 eV. 



solutions have been described in 32]. Our intent in this 
section is to investigate, via a nonlinear analysis in phys- 
ical space, the formation and evolution of neutrino lumps 
and their properties in redshift, in order to estimate the 
final gravitational potential characterizing the lumps as 
a function of their final scale. 

With this aim in mind, we solve the nonlinear Navier- 
Stokes equations in an expanding Universe and position 
space with comoving coordinates x: 

5l = -v v -V5 v -(l + 5 1/ )V-v v , (13) 
v' u = -(H-04') v v -(vu-V)v„ 

+V($„ + 0ty) , (14) 

A5c/> = -pa 2 6„p„ , (15) 

^ v = ~5 v p v . (16) 

Here p v is the background neutrino energy density and 
5 U = Sp u /p u is the relative neutrino density perturbation 
(~ 1 when reaching nonlinearity) . The vector v„ is the 
velocity for neutrinos. More precisely, it describes the 
peculiar comoving velocities - it vanishes for neutrinos 
with constant comoving coordinates. The evolution of 
the velocities is driven by the gradients of the gravita- 
tional potential and cosmon field, with the usual Hubble 
damping and quadratic term arising from particle num- 



ber conservation. The velocity dependent term /3</>'v„ 
in ea. (IT4|) is not present in the standard Navier Stokes 
equations. It accounts for momentum conservation, re- 
flecting the fact that the neutrino mass changes in time 
as m'/m = — 04' ■ This friction term can be rigorously 
derived within the fully relativistic equations, as outlined 
in the last section of this paper. 

Equation (fl~6)) is the Poisson equation for the gravi- 
tational potential that we have indexed as to clarify 
that it only comprises the neutrino contribution. (The 
sign convention for the gravitational potential in eqs. (|14[) 
(ITBl matches the linear equations in 20] with A corre- 
sponding to — k 2 in momentum space). Eq. (|15p is the 
perturbed Klein Gordon equation in the limit in which 
time derivatives are negligible with respect to the spatial 
ones, as it holds in the Newtonian approximation. Effec- 
tively, eq. (|13p relates the time evolution of the neutrino 
ovcrdensity 8 V to the divergence of its corresponding mo- 
mentum density. The time dependence of the momen- 
tum density itself is directly connected, via eq. (fT4|) . to 
the given forces: in addition to the gravitational force 
F g = V$„, the cosmon mediated fifth force F = f3\7S(f) is 
present, in the form derived already at a linear perturba- 
tion level (see sec. IVII[) . Combining eqs. (fT5"|) and (TIT))) im- 
mediately shows that 84 ~ 2/3$, clarifying that the scalar 
field mediates a force of order |F| = \0VS4\ ~ 20 2 \F g \. 
For a choice of ~ —50 this is about 5000 times stronger 
than gravity. 

In section IVlIl we will show that eqs. (IT3|) and (fl"4f can 

be obtained by considering the appropriate limits of the 
fully relativistic equations derived from the Bianchi iden- 
tity in presence of an external source [35j 

V 7 T; = Q M = ~T^<t> , (17) 

where TJJ is the stress energy tensor of the neutrino fluid. 

For convenience we display eqs. ([T3"l - 1TB ]) also in terms 
of cosmic time t and physical (not comoving) coordinates: 

dp 

— f = -V(p i/ V to t) - 04 Pv - PV(j>p u Vtot; (18) 

ot 

^ = (Vvtot + 13 4>) vtot + V ($„ + 04) (19) 

- (vtofV) Vtot , 

Here p v and 4 are the local neutrino density and cos- 
mon fluid including the fluctuations and v to t is the total 
neutrino velocity, composed of Hubble flow and peculiar 
velocity as 

dv 

vtot = -j- = Hr + v, . (20) 
at 

Equations (fT5]l and (jTSJ) may be combined to yield the 
conservation equation of the momentum density p = 
p u vtot, namely 

Pi = fdisp(i) + fattr(j) , (21) 
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where we have defined 



fattr(i) = Pudi(^ v +j3(/>) , 



(22) 



The attractive force f a tt r may eventually be balanced by 
the countering force associated to the velocity dispersion, 
fdisp- We have omitted in eq.(|22|) and eq. (fl~8|) an addi- 
tional pressure force 



fprcss(i) = -3(i? - P<f)p u ^tot(i) 

Pu = w v {p v + 8p v ) , 



(23) 



which is important only when neutrinos are still relativis- 
tic. It prevents the growth of neutrino overdensities for 
z > 10. 



IV. EVOLUTION OF NEUTRINO LUMPS 
A. Method 

We integrate the set of nonlinear equations (I13H16[) on 
a 128x128x128 point 3-dimensional spatial grid of fixed 
size L 3 , by means of a method of lines with an adap- 
tive time stepper. The length L is chosen to match the 
given initial profile, such that a good resolution around 
its maximum can be obtained. We impose cubic periodic 
boundary conditions, p(x) = p(x + L) with Li = 0, L. 
For each time step, the Poisson equations (TT51) (TT61 are 
solved employing a Fast Fourier Transform routine. 

We start by considering the formation of a single lump 
within our box. As an initial density profile we consider 
a gaussian: 



6 u (x) = hi 



2 / 2 



(24) 



The initial density amplitude hi n is chosen to match the 
corresponding linear overdensities at the matching red- 
shift z m , an issue which will be discussed in more detail 
in the next section. A convenient width of the box cor- 
responds to about L ~ 6r in . We start the numerical 
integration early enough (z ~ 9) in order to allow us 
to start with rather arbitrary small neutrino velocities 
Vi n (see below). Independently of their precise values 
there is enough time for their radial component to adapt 
such that linear perturbations are matched by the time 
we reach the range of validity of the nonlinear equations 
specified above. We also investigate alternative initial 
conditions for the velocities. The dependence of our re- 
sults on initial conditions is discussed in the next para- 
graph, while illustrating our results. 



first phase in which the lump expands with the back- 
ground, peculiar velocities will take over. The lump will 
then decouple from the background expansion and con- 
tract with increasing inward velocities. A useful quantity 
to illustrate the rapid increase in velocity is the kinetic 
energy of the lump, defined in physical coordinates as 



El 



&Pv v t 2 ot d 3 a 



(25) 



with total velocity v to t = v„ + "Hx = v„ + Hr, being 
the sum of peculiar and expansion contributions. Com- 
paring the kinetic contribution to the potential energy 
can give us an indication of when and how fast the lump 
approaches virialization. The potential energy is defined 
as 



Epot = -\& J 5p y 5(j)d 3 x 



(26) 



where V denotes the physical volume of the overdensity. 
We have omitted here the small gravitational potential 
energy. Both E^in and E pot are evaluated in practice 
in cartesian coordinates rather than in spherical ones, 
as late phases of the collapse might involve important 
anisotropies. 
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FIG. 2. Linear (solid, red) and nonlinear (long-dashed, green) 
neutrino overdensity <5„ vs redshift, until virialization occurs. 
The nonlinear overdensity is evaluated in the center of the 
lump. The ratio of kinetic over potential energy associated to 
the lump is also shown (short-dashed, blue). The comoving 
(physical) initial lump radius r^ n = 45 Mpc (4.46 Mpc) fixes 
the box size L — 270 Mpc for the simulation and corresponds 
roughly to a scale k/h — O.lMpc -1 for h = 0.72. It also 
determines a final physical radius of the lump Rf = 3.21ri„ = 
14.3 Mpc. 



B. Turn over and free fall 

The nonlinear evolution of a perturbation as resulting 
One of our aims is to understand and illustrate the evo- from the integration described above is shown in Figf2] 
lution of neutrino lumps in time: we expect that after a It is compared to the corresponding linear evolution. We 
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FIG. 3. Evolution of potential (pink, dotted) and kinetic en- 
ergy (red, solid) versus redshift. The two main contributions 
to the kinetic energy are also shown: the energy associated 
with the Hubble expansion is depicted in green (long-dashed) 
and the one due to peculiar motion in blue (short-dashed). 
The size of the lump is the same as in fig[U An energy of 
10 erg corresponds to roug hly 6 • 10 77 eV. 



have chosen n n = 45 Mpc (4.46 Mpc) in comoving (phys- 
ical) units and find a characteristic final size of the lump 
Rf = 3.21rj n = 14.3 Mpc in physical units. The precise 
definition of Rf will be given below. The ratio of kinetic 
to potential energy is also shown in the plot. The evo- 
lution of the potential energy and of the kinetic energy, 
together with the split of the kinetic energy into expan- 
sion and peculiar components, is further shown in Figj3] 
for the same final scale. For the quantitative evaluation 
of the energies we choose a volume that extends to a ra- 
dius where the density contrast reaches (1/50) times the 
central density contrast. 

By looking at Figj2] and Fig|3j we can identify three 
redshift ranges: 

: [i] 1 + z > 3: In this range, density perturbations 
are still linear; total velocities and therefore ki- 
netic energy are dominated by the expansion term 
8p l/ V?-x 2 (= 8p v H 2 r 2 ) present in eq. flU} while the 
peculiar velocities v are negligible. At this stage, 
both potential and kinetic energy are increasing due 
to the accretion of more neutrinos in the lump and 
to the increase in the neutrino mass. The potential 
energy E pot increases somewhat faster (both 5p v 
and 8(f) are proportional to m„), which results in a 
decreasing ratio E kin /E pot . 

: [ii] 2. 5 < 1 + z < 3: At z ~ 2 the attractive fifth 
force starts to become the dominant contribution 
in eq. (|2ip . In this regime we can trust the non- 
linear evolution for the chosen scale. The nonlin- 
ear density perturbation detaches from the linear 



one. Also the radial peculiar velocities start to be- 
come nonnegligible, adding an inward contribution 
to the outward expansion term. As a consequence, 
for 1 + z < 3 the increase of the kinetic energy is 
slowed down. The latter effect is more pronounced 
for inner shells of the profile, for which the radial 
velocities change sign, while outer shells still tend 
to expand with the background Universe. During 
this period E pot keeps increasing, with a slightly 
steeper slope than in the previous range due to the 
effects of nonlinearity on the neutrino density per- 
turbation. 

: [iii] 2.3 < 1 + z < 2.5: The nonlinear effects on the 
density perturbation become very pronounced. Pe- 
culiar velocities dominate over the expansion term 
and induce a rapid increase of the kinetic energy, 
due to the inward velocities. As a consequence, the 
characteristic scale of the lump shrinks until virial- 
ization is reached. Equation (|14p is now completely 
dominated by the fifth force. The lump has effec- 
tively decoupled from the background and evolves 
according to the laws of free fall. At the end of this 
period the tangential peculiar velocities as well as 
irregularities in the spherical symmetry of the lump 
grow rapidly. The borders of numerical precision 
are reached. The latter effects lead to a highly 
anisotropic behavior, gradients become very large 
and velocities change directions almost randomly. 
The timestepper eventually stops the integration. 
Note that our code stops as the ratio \Eki n /E pot \ 
approaches the value of 1/2, corresponding to the 
virialization condition. The latter has not been in- 
troduced by hand in the equations but is automat- 
ically reached due to increasing tangential peculiar 
motions and deviations from spherical symmetry. 

As the neutrino lump undergoes the different evolu- 
tionary stages, the density profile adapts according to our 
equations once the nonlinear regime is reached. While 
density perturbations are still linear, the profile expands 
with the background, leaving its shape almost inaltered. 
The only major change is the buildup of a surrounding 
underdensity as neutrinos flow into the lump. Note that 
for a single lump the underdensities are always much 
smaller than the background density, even when the over- 
densities become large. When collapse is approached, the 
profile changes its shape. As outer shells turn around at 
smaller redshifts, the slope of the profile changes. Its 
initial gaussian shape changes into a profile which resem- 
bles common density profiles for dark matter halos, such 
as the Navarro-Frenk- White profile. We plot the density 
profile for different redshifts in figj4] 

The size of the neutrino lump is the most characteris- 
tic quantity which distinguishes the evolution of different 
types of inhomogeneities. If we start with vanishing ini- 
tial peculiar velocities and adjust hi n (rj n ) according to 
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FIG. 4. Neutrino overdensity profiles at redshifts z = 8.9 (red, 
solid), z = 2.2 (green, long-dashed) and z = 1.3 (blue, short- 
dashed) in terms of physical distance r from center, averaged 
over all angles. The central overdensities are So = 5(0, z) — 
1.1-10" 5 , 1.6-10 -2 , 250, respectively. To illustrate the change 
in shape we have normalized the profiles by dividing out So. 
The radius is given in physical units. An initial broadening 
due to the expansion is followed by a strong concentration 
process. An apparent discrepancy between the width of the 
profiles and the radii shown in fig[5]is caused by the different 
definition (f27)) . 
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FIG. 5. Radii of the neutrino overdensity profile at different 
amplitudes, for a final maximal radius of Rf — 14.3 Mpc. We 
plot the time evolution of the physical radii where the neu- 
trino overdensity reaches a certain fraction of the core density. 



the linear evolution of small perturbations, our initial 
conditions depend on only one relevant parameter ri n . 
We are interested to determine the final size Rf of the 
neutrino lump in terms of r;„ . An understanding of the 
time evolution of the different length scales in the lump 



FIG. 6. Late evolution of different radii for a final maximal 
radius of Rf = 14.3 Mpc. This plot zooms fig. [5]to a smaller 
range of z, such that the more irregular behavior close to 
virialization becomes visible. 



is also important for the precise definition of the volume 
integral in eqs.([25]) (|26|) . 

It is therefore interesting to consider different radii 
of the profile, characterizing shells at different distance 
from the center of the lump. In order to define radii 
corresponding to inner or outer shells, we consider the 
physical radius i?(A, z) for which the amplitude reaches a 
given fraction A of the central amplitude: 5{z, R(X, z)) k, 
X5(z,0). The radius R(X, z) depends on the redshift z. 
Higher values of A correspond to amplitudes which are 
closer to the central value and therefore to inner shells 
inside the lump. Since the lump has no spherical symme- 
try we have to make the definition of R(X, z) more pre- 
cise by taking the largest value of r for which the density 
A<5(0) is reached, or formally define the radius R(X, z) as 
the superior limit of the {r} ensemble. In formulae we 
have: 

R(X,z):=sup{\r\ \ 8 u {z,r) > A«S„(z,0)} . (27) 

Here the factor A is a constant < A < 1 and all spatial 
coordinates f are again physical. 

We show the time evolution of R(X, z) in units of rj„ 
in FigsE] and HI We have chosen four different values 
for A, i.e. A G {0,1/50,1/4,3/4} and display R (z) = 
R{\ = 0, z), R 1/50 (z) = R{^,z), Ry A {z) = R(\,z) and 
i? 3 / 4 (z) = R(j, z). As expected, for larger values of A the 
corresponding shell turns around earlier. For compari- 
son, we have also plotted the slope corresponding to the 
scale factor a{z). As expected it corresponds to the slope 
of R\ within the redshift regime in which the expansion 
is still dominant. As a measure of the size of the lump 
we employ Ri(z) = Ri/ 50 (z). We assume that the change 
of R after virialization can be neglected, such that the 
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final size Rf is given by Rf = R{z V i r ) = i?(l/50, 

C. Onset of virialization 

The late evolution of the different radii displayed in 
Fig|6] shows a rather irregular behavior. This indicates 
the onset of virialization, as also suggested by the ra- 
tio Ekin/Epot in Fig. ^j. The reason is that transversal 
velocities become important in this redshift range. In 
order to illustrate this we analyze the evolution of non- 
radial velocities within the lump. In fig. [7J we plot the 

ratio -Ekin, nonradial/ -E"kin, peculiar , with 

Skin, nonradial = ^ J $pu 5^1 d 3 X , (28) 
#kin,pcculiar = 7. [ ^Pu^l d 3 X , (29) 

^ Jv 

where <5v = v — (R' /R)r, v is the peculiar velocity, R the 
radius of the lump and r is the radial coordinate. While 
for higher redshifts the nonradial contribution to the pe- 
culiar kinetic energy is entirely negligible, it increases 
rapidly once the evolution of the lump has predominantly 
turned around. 

In order to understand this, it is useful to consider 
a shell with a given radius R s (t) and analyze equa- 
tion (|13l) to see how deviations from its mean velocity 
evolve 36]. For this purpose, we transform into a coor- 
dinate system in which the mean of the shell is at rest, 
x = R s (t)t = (a/i?)x, where r and x are the phys- 
ical and comoving spatial coordinates, respectively. We 
then decompose the peculiar velocities into a radial part 
V|| = x' and a nonradial part <5v s and insert it into equa- 
tion (|14[) . Appropriately transforming spatial derivatives 
leads to the expression. 

K = - (§ - PA *v s (30) 

-TT [(Sv s -V)Sv s - V_L($„ + /3<ty)] • 

The gradient Vj_ in front of the potentials shows that 
only nonradial components influence the evolution. This 
may also be understood by noting that the underlying 
situation is formally equivalent to the evolution of the 
velocity perturbation within a Universe expanding with 
the rate R'/R. The structure should thus simply resem- 
ble equation ([T4")) with an altered Hubble function. The 
additional factors of a/R are due to the choice of time 
coordinate. 

As long as the shell is expanding, the friction term 
R'/R behaves as a damping force. Once the shell turns 
around, however, R' becomes negative, and the term ef- 
fectively enhances the build-up of nonradial velocities. 
Due to the large modulus of the collapse rate, this hap- 
pens quite rapidly. Of course, an exactly radial flow will 
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FIG. 7. Ratio of the nonradial part of the peculiar kinetic 
energy and peculiar kinetic energy as a function of redshift. 



remain radial by virtue of rotation symmetry (up to nu- 
merical errors). We have therefore started initially with 
small anisotropies in the velocity (or density) distribu- 
tion. For this purpose, we have imposed an additional 
initial velocity field v ran d, whose directions were ran- 
domly distributed. Its amplitude at each space point 
|v ran( j| was chosen to be of an order of a few percent 
of the mean initial expansion velocities. We find that 
the late stages of the flow are independent of the precise 
choice of initial conditions for the nonradial velocities. 



V. GRAVITATIONAL POTENTIAL OF 
NEUTRINO LUMPS 

A. Single lump potential 

The key quantity of interest is the characteristic gravi- 
tational potential of a neutrino lump. This will influence 
the peculiar velocities of galaxies or gas, or the CMB- 
anisotropies via the ISW effect. Of course, the gravita- 
tional potential depends on the distance from the center 
of the lump as shown in figlU We have indicated the scale 
Rf = Ri/ 50 (z vir ) by a dot. FigfS] demonstrates that the 
radius R1/50 i s a typical scale which characterizes the 
gravitational potential of the lump. We may define 
as the gravitational potential in the center of the lump, 
^ V: R t as its value at the radius Ri defined in eq. (|2"7|) . and 
^v.Ri/3 a corresponding value at distance Ri/3 from the 
center. We plot in Fig |H] the time evolution of the gravi- 
tational potential, for a lump with final radius Rf = 14.3 
Mpc. Note again that & u is obtained as the solution to 
the Poisson equation (fl"6j) . 

In fig 1 101 we consider families of lumps with different 
initial conditions and show the gravitational potential of 
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FIG. 8. Dependence of gravitational potential on distance 
from the center at redshift of virialization z = 1.3. The dot 
indicates Rf. 



FIG. 10. Gravitational potential as a function of the final 
size R of the neutrino lump (at redshift z = 0). We define 
7? as the radius at virialization Rf, assuming that it remains 
approximately constant from virialization up to now. 
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FIG. 9. Neutrino gravitational potential as a function of 
redshift for a fixed final scale Rf = 14.3 Mpc. We evaluate 
$„ for different distances from the center of the lump, namely 
at Ri = -R1/50, Ri/3 and in the center of the lump ($i/, c ). 




ally merge to bigger structures, such a merging process 
can be very different from the essentially spherical in- 
fall investigated in this paper. One typically expects a 
smaller gravitational potential for a big structure aris- 
ing from merging as opposed to spherical infall, since the 
"substructures" carry tangential velocities. In this re- 
spect our result for <£>„ should be considered as an upper 
bound, in particular for large values of R. It should be- 
come a good approximation for small enough R where the 
possible merging processes do not play a dominant role. 
Without extended numerical simulations it is difficult to 
assess the value of R beyond which the true character- 
istic gravitational potential remains substantially below 
the values in figHOl With these words of caution the final 
value for a virialized structure of typical size ~ 10(100) 
Mpc is - 10~ 6 (10- 4 ). 



B. Neutrino induced cosmological gravitational 
potential 



the lump versus its final size Rf. We plot the nonlinear 
value of the gravitational potential $„ as obtained from 
the numerical solution of eqs. (|13fll6]l for the definition of 
Rf = R(l/50, z V i r ) according to eq. (f27j) . This value has 
to be handled with care. We have investigated in this 
paper only the situation where a single lump forms, es- 
sentially independently of possible surrounding lumps. In 
our numerical work, this follows from the choice of initial 
conditions. In the true Universe, the situation may be 
much more complicated. Many lumps of different sizes 
are expected to form, competing for the available neu- 
trinos. Even though these "initial lumps" may eventu- 



For an estimate of the cosmological impact of the grav- 
itational potential of the neutrino lumps one needs a rela- 
tion between the characteristic potential for a single lump 
and the average cosmological value of the gravitational 
potential at a given scale, expressed as the Fourier com- 

(c) 

poncnt I . One expects a typical suppression factor 7 C 
for the cosmological gravitational potential as compared 
to the potential of a single lump with size Rf = Tr/k, 

l$ (C ll 7T 

^ s ^n-* = s ■ (31 » 

I v,R\ 
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For clarity we here indicate the cosmological potential 
and the single lump potential by superscripts (c) and (/). 
If all lumps would have the same size, the suppression 
should be roughly given by the fraction of the volume 
occupied by neutrino lumps. For a rough numerical esti- 
mate of 7 C we may then use the relations 

_ ^y(lumps) 



t J y v lumps 



pl c) V hoi = N v 



(32) 



with Vi unl p the total volume occupied by neutrino lumps 
and Vhor the volume within the cosmological hori- 
zon, V hor « ^ (3000 Mpc) 3 , and define by = 
jyOumps)^^ ^ e total fraction of neutrinos which arc 
bound in lumps of size R = n/k. This estimates 7 C in 
terms of the ratio of the average energy density pi of the 
neutrinos in the lumps and the cosmological background 
density pi 



Vi 



umps 



v ho 



Pu 



(33) 



For a typical neutrino density averaged over the volume 
of the lump within radius R = Ri / 50 we find pi /pu ~ 
100. In our simulation we find typical values for Fu 
for single lumps in a box to be « 0.1 at 2„irj but 
this depends of course on the size of the box. Even if all 
neutrinos are finally bound in some lump one typically 
has F$r (fc) < 1 since not all are in lumps of size n/k. We 
will take F$> 1) = 1/4. 

For an alternative rough estimate of "f c we have ran- 
domly distributed virialized neutrino lumps with a fi- 
nal profile according to our numerical solution over a 
cosmological volume V c . The number of lumps was 
chosen such that the total number of bound neutrinos 
matches F„ N v . The corresponding gravitational poten- 
tial <&£f^(x) has then been Fourier transformed in order 
to extract <by k = ^- J d 3 x $^(x) exp (—• ikx). We show 
in fig JTT] oversimplified distributions where all lumps cor- 
respond to a single size R. Nonetheless, using eq. (I3T1) . 
the order of magnitude of 7 C can be extracted from figOU 
We show the resulting 7 c (fc) in figH2l together with the 
estimate (f33|) . for Fu = 1/4. We also show a mixed 
distribution of lumps with final radii of Rf = 14.3 Mpc, 
31.8 Mpc and 63.6 Mpc, with equal numbers of neutrinos 
in each sort of lumps and a total F„ = 1/4, as well as 
an analogous mixed distribution where the mass of the 
lumps is concentrated in the center (point masses). 

A suitable interpolation of figfT2l as indicated by the 
dot-dashed light blue line, corresponds to a simple fitting 
formula 



<fi c) (fc) 
fc = 3.9 



1Q- 2 h/Mpc,K = 2 



(34) 



This fit should only be taken as a useful order of mag- 
nitude estimate. FigfT2l illustrates that the final order of 



magnitude of the gravitational potential is quite sensitive 
to the choice of distribution. The expected cosmologi- 
cal distribution of neutrino lumps involves a distribution 
over a substantial range of lump sizes. Moreover, large 
lumps may have substructures. This makes it hard to es- 
timate the relevant number Nu m (k) without explicit 
cosmological simulations of structure formation. On the 
one side, the total fraction of neutrinos found in lumps 
may come close to one. On the other hand, the distri- 
bution of the available neutrinos over lumps with various 
sizes may reduce the effective fraction Fv relevant for 
a given scale. Especially for small k near the horizon 
the clumping may not have had enough time to bind a 
large fraction of neutrinos. Small F^ (fc) would result in 
a further suppression of the effective 7 c (fc). Furthermore, 
note that our fit is likely to overestimate the power on 
small scales, i.e. for large fc. At these scales, the exact 
shape of the profile is still of substantial importance, as 
can be seen from fig |12l However, the simple fitting for- 
mula (|34|) with k = 2 is a good estimate at large scales 
and therefore a viable mean to illustrate the cosmological 
implications of neutrino lumps. 
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FIG. 11. Gravitational potential as obtained by Fourier trans- 
forming distributions of lumps of given size. For comparison 
we have included the single lump gravitational potential at 
distance R — Ri/50 from the center. 



In figiT3] we plot |$^| as a function of fc according 
to eq. (IMl) . The substantial uncertainties in this esti- 
mate are reflected by error bounds which correspond to 
a multiplication of 7 c (fc) by factors of 4, corresponding to 
Fu — 1, and 0.1. The upper limit may be regarded as 
a rather solid bound - effects of merging dynamics and 
Fu < 1 reduce <&„(£;) as compared to this bound. A 
more reliable estimate of $ u k needs the actual distribu- 
tion of neutrino lumps as a function of size and mass - 
this may perhaps require extended A^-body simulations. 
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FIG. 12. Suppression factor y c {k) as obtained by using the 
definition (|3ip for distributions of lumps of size Rf = 14.3 
Mpc (solid, red), Rf = 63.6 Mpc (long-dashed, green) and a 
mixed distribution (short-dashed, blue) . The dotted pink line 
marks the gravitational potential of the mixed distribution of 
point masses, while the dot-dashed light blue line corresponds 
to a suitable interpolation. The estimate (|33p is also shown 
(double-dashed, black). It almost coincides with the interpo- 
lation of the mixed point mass distribution. The error bounds 
displayed (shaded, light blue) correspond to 47 c and 0.1y c . 



In fig|T!|]we have also included the strongest possible up- 
per bound for the cosmological gravitational potential, 
obtained by assuming all neutrinos in the horizon are 
bound within a single, point-like lump. This results in 
the following expression: 



2fc 2 " yu '2M 2 k 2 
m v {t Q ) 1.07- 10- 8 (h/Mpc) 2 
2 eV P 



(35) 



Any physically realistic scenario will lead to substantially 
lower gravitational potentials. 

For scales R > R PN ~ 500 Mpc, k/h < 6.3 • 10~ 3 
Mpc , the Newtonian approximation breaks down and 
our estimate of the nonlinear <5> v would have to be cor- 
rected by non-Newtonian effects. The exact matching 
of the linear evolution of $„ for scales outside the hori- 
zon and the nonlinear results would require a nonlin- 
ear analysis which also includes post-Newtonian terms. 
This unknown region 10 -3 ^ k/h ;$ 6 ■ 10 -3 may leave 

the strongest observational imprint by the ISW effect. 

FigJI51 can be taken as a clear illustration that the for- 
mation of nonlinear structures makes the output of lin- 
ear Boltzmann codes unreliable when density fluctuations 
become nonlinear. Instead of increasing continuosly, the 
gravitational potential will saturate at the typical val- 
ues found by the nonlinear analysis. An unjustified ex- 
trapolation of the linear approximation overestimates the 



FIG. 13. Gravitational potential as a function of the scale k 
obtained from the nonlinear (red, solid) evolution. The er- 
ror bounds (shaded, red) arise from the bounds for f c . For 
scales k smaller than fcpjv (solid square) the Newtonian ap- 
proximation is no longer valid and the extrapolated nonlinear 
potential is drawn as a dotted line. The dot-dashed black line 
corresponds to the extreme bound (1351) . 



gravitational potential by many orders of magnitude for 
scales of a few hundred Mpc. For the parameters used 
in this paper, m v — 2.1 eV, /3 = —52, the linear extrap- 
olation would predict enormous oscillations in the spec- 
trum of CMB anisotropics for angular momenta I ~ 100, 
very strongly in contradiction with observation. Placing 
reasonable bounds on the gravitational potential, consis- 
tent with our findings for the nonlinear evolution, com- 
pletely eliminates this feature and makes the spectrum 
of anisotropies insensitive to growing neutrinos in this 
range of I. 

On the other hand, in the range of smaller I ss 10, 
a large time variation of the gravitational potential for 
scales 10~ 3 h/Mpc < k iS 10~ 2 h/Mpc could leave a big 
imprint on the CMB spectrum. A too large ISW ef- 
fect could rule out such a scenario. We conclude that in 
principle the ISW effect provides a viable mean to con- 
strain growing neutrino cosmologies, at least in the case 
of large neutrino masses. For smaller neutrino masses 
the neutrino fraction Q u of the energy density is smaller, 
thereby presumably reducing the ISW signal. 



VI. INITIAL CONDITIONS 
A. Matching of the radius 

The size of a given neutrino lump under investigation 
is determined by the parameter fj„. For a comparison 
with linear perturbation theory, as in fig llOl we have to 
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relate r;„ to an appropriate comoving wave number k. 
We will see in the next section that for a given scale 
there is a range of redshifts at which the fully relativistic 
linear and the Newtonian nonlinear evolution coincide. 
We can make use of this in order to approximately fix 
the relation between our initial width parameter ri n and 
the comoving momentum scale k for the associated linear 
fluctuation. For this purpose we choose a redshift z m 
within the discussed range. At the given redshift z m we 
choose ri n such that the nonlinear and linear quantities 
relate via 



7T 
k 



(36) 
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Note that the exact choice of z m is not important. Within 
the linear regime R/a — const , since the radius scales 
with the background scale factor. Hence, different match- 
ing redshifts yield equivalent results for the relation be- 
tween Ti n and k. 



FIG. 14. Dependence of the gravitational potential at virial- 
ization and virialization redshift Zvir on the initial amplitude 
hin for ri n — 45 Mpc. The oscillations in 3>„ reflect the oscil- 
lation of the neutrino mass as a function of z V i r . 



B. Matching of the amplitude 

Taking the radius of the overdensity as a free param- 
eter, the amplitude hi n still needs to be fixed. From 
a primordial gaussian fluctuation spectrum we expect a 
characteristic distribution of initial amplitudes around 
zero, with a mean square dispersion given by o~u n (ri n ). 
Actually, we find that the redshift of virialization z V i r de- 
pends monotonically on /i in , while the final gravitational 
potential of the lump remains rather insensitive with re- 
spect to the precise choice of hi n . We demonstrate this 
in figH3J This approximate 'universality' of the lump 
properties after collapse may be of substantial help for 
observational investigations of our scenario. 

A good estimate for a characteristic amplitude hi n can 
be obtained as follows. We consider the overall perturba- 
tion <J {fin, hi n ) defined as the convolution of our initial 
gaussian density configuration <5(x) and a window func- 
tion of size r in : 

cr 2 (r in ,h in ) = -^ [ \6(k)\ 2 \W(k)\ 2 d 3 k , (37) 
v Jv 

where <5(k) and W(k) are the Fourier transforms of the 
density configuration <5(x) and of the window function 
W(ri n ) respectively. The best estimate for hi n corre- 
sponds to the value for which a(ri n ,hi n ) = aun(rin)- 
Here au n (ri n ) is computed by using in eq.(|37|) the same 
window function W(ri n ), but now the linear fluctua- 
tion spectrum 5u n {k) instead of the Fourier transform 
of eq. ([24|) . The linear Su n (k) are evaluated from a 
Boltzmann code for growing neutrinos and taken at red- 
shift Zi n . Note that the matching is done in fc-space 
and therefore requires no explicit calculation of the re- 
lation between length scales r and momentum k. Since 



c 2 (^m, hi n ) ~ hf n we obtain a simple estimate for a char- 
acteristic initial amplitude. The dependence of the ini- 
tial amplitude h in on the initial radius r in is displayed in 
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FIG. 15. Relation between initial radial parameter r% n in 
comoving units and dominant initial amplitude hi n . 



Our choice of hi n neglects two effects: the dragging 
of neutrinos by dark matter fluctuations and the nonzero 
pressure at early stages of the evolution. We have checked 
that these effects largely cancel. This is demonstrated by 
the good agreement with full linear perturbation theory 
(including both effects) in the matching range, as dis- 
cussed in the next section. 

In concluding this section, we mention a series of checks 
that have been done in order to verify the stability of our 
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results. 



• We have checked how the formation of the lump 
depends on different initial velocities. In particu- 
lar, increasing initial nonradial velocities leads to a 
faster collapse. All features, such as size and total 
mass of the final lump remain unaltered. In partic- 
ular, the final gravitational potential is unaffected 
by a change of initial velocities. 

• We have investigated the formation of two lumps 
within our box, both in the case in which they are 
far enough apart to evolve independently and in 
the case in which they are close enough to merge. 
For a given final size Rf of the lump, neutrinos 
distribute equivalently in the two cases, with no 
relevant effect on the gravitational potential. This 
scenario does not cover, however, the interesting 
and probably realistic case of many lumps moving 
for a while under their mutually attractive cosmon 
forces. 

• We have considered the case of an initially non- 
spherical profile, with a pronounced ellipticity as 
in the following expression 



For fast comparison to the nonlinear equations we re- 
call here the set of linear perturbation equations specified 
in the case in which the Newtonian limit applies. (For 
the full, closed, system of linear perturbation equations 
see |20| , where linear perturbations for growing neutrinos 
have been extensively investigated.) They read 



8'^ = 'H-L(w ( f > - 0^)6^ - (1 + w<f,)kv4 

[(l-3<)<5 -(l-3cl)6 v ] , 

P<f> 



(38) 



S' V = 3(H- P<f>') {w v - cl)8 v - (1 + w v )kvv • (39) 

In deriving these equations, we have defined the line 
element as given by ds 2 — — a 2 [(l + 2^)dr 2 — (1 + 
2$)7yda; l dflf J '], where $ corresponds to the usual grav- 
itational potential and we are working in Newto- 
nian gauge. The equations for the density contrasts 
5i(k) = y J 8i(x)exp(—i\t ■ x)d 3 x (defined as the Fourier 
transform of the local density perturbation <5j(x) = 
5pi(x)/pi(x) over a volume V) involve the velocity per- 
turbations, which evolve according to 

vL = -H(l - Sw^v^, - (3(f>'(l - Zw v )—v^, 



8 v {yi) = K M e w ™ *"* , 

with w z S> w xy . Again, for a given Rf the results 
do not seem to be affected in an important way. 

• Finally, we have estimated the effect of Fermi pres- 
sure within a Thomas-Fermi approximation as a 
function of the lump size. The ratio of Fermi force 
and fifth force, increasing with the scale, never hap- 
pens to be bigger than a few percent and there- 
fore has a negligible effect on the estimated gravi- 
tational potential. 



VII. RELATIVISTIC EQUATIONS 

In this section we further discuss the range of valid- 
ity of the numerical integration and clarify the choice 
of initial conditions. We are interested in the range of 
scales and redshifts where neutrino perturbations are ex- 
pected to grow nonlinear as found in [20j |. when neutri- 
nos are nonrelativistic and pressure terms are negligible. 
The integration of the nonlinear equations starts inside 
the linear regime, when neutrinos start becoming non- 
relativistic, at Zi n ~ 9. At this time the system is still 
correctly described by the set of fully relativistic linear 
perturbation equations, in which neutrino pressure terms 
are included. We want to make sure to have a proper 
matching of the regime where we can trust the nonlinear 
equations used in the numerical solution to the regime 
where the relativistic linear equations are valid. 



1 + Wcj, n + 

r , , kn T + Pk8(j> — , (40) 

3 1 + W0 P<h 1 



P0 1 + W<f, 



v' v = (l-3«v)^-%- 



1 + w v 



kct 



fik8<j> 



1 + w v 
1 — 3uv 



, , 2 w v 

fc * _ q TT TTf 

3 1 + w v 



1 + W v 



(41) 



These equations reduce to the linearized Navier Stokes 
equations if we set neutrino pressure terms to zero and 
consider the case of no anisotropic stress for which = 
— thus obtaining for the neutrino component: 



5' 



= —kv v 

l = -(n 



(42) 
(43) 



The latter equations, here written in Fourier space, are 
equal to the ones obtained by linearizing the Navier 
Stokes equations (jlUl ITU) . 

For a certain first time period after we start the in- 
tegration of the nonlinear equations, neutrino pressure 
terms will not yet be negligible and therefore the output 
of the linear equations (|39I41[) will differ from the out- 
put of the nonlinear equations (|13H16[) . in which we do 
not include pressure terms. We will therefore consider 
the nonlinear equations without pressure terms as reli- 
able only after a certain redshift z n i, while prior to zi we 
trust the relativistic linear equations. There is a redshift 
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range where both equations are valid. We match the evo- 
lution at a certain redshift z m in the range z\ < z m < z n \ , 
which depends on the size of the lump under investiga- 
tion. This matching is done by choosing appropriate ini- 
tial conditions for the numerical solution. More precisely, 
we require for the 'matching' redshift z m {k) that pressure 
terms effectively become negligible. For z < z rn we will 
then consider the output of the nonlinear integration as 
valid. 

In fig [in] we show the critical redshift z cr for which 
neutrino fluctuations become well approximated by pres- 
sureless nonrelativistic particles. The critical redshift 
is plotted as a function of the length scale defined as 
R/a = 7r/fcj where k is the momentum scale. We define 
the critical redshift as the value of z at which the rela- 
tive size of pressure terms drops below a certain value in 
the linear equation for v v ([41~j). The shaded region cor- 
responds to relative pressure contributions between 1% 
and 10%. Above the upper (green) line pressure contri- 
butions are > 10%; below the bottom (blue) line pressure 
contributes for less than 1%. 
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FIG. 16. Critical redshift versus momentum scale. The crit- 
ical redshift is defined as the z at which pressure terms be- 
come negligible in the equation for v v , for each scale. The 
blue (dark shaded) region corresponds to relative pressure 
contributions less than 1%; the grey (intermediate shaded) 
region corresponds to relative pressure contributions between 
1% and 10%. The above light grey region corresponds to 
pressure contributions > 10%. 



Our nonlinear equations are meaningful for those scales 
and redshifts at which radiation, baryons and CDM den- 
sity perturbations are subdominant with respect to neu- 
trino density perturbations, as it happens for large scales 



at late times due to the rapid effect of the fifth force act- 
ing on neutrinos only. Only in this regime it is justified 
to restrict the contributions to the gravitational potential 
in (TT6"|) to neutrinos. Finally, note that although lumps 
are expected to form at very large scales (~ 200 Mpc), 
these are still well within the horizon, which justifies a 
Newtonian approach. 

In closing this section, we would like to mention that 
it is also possible to show that the nonlinear equations 
described by ((T3l - [l~6l) follow from the set of full relativis- 
tic equations taken in the Newtonian limit. We do not 
intend to show this explicitly here, but just outline the 
basic procedure. For the derivation of the nonlinear rela- 
tivistic equations we consider the case in which there is no 
anisotropic stress, *5> = — <& so that the line element takes 
the form ds 2 = -a 2 [(l ~2$)dr 2 -(l+2^)j ij dx i dx j ]. We 
consider the equations in the Newtonian limit, restrict- 
ing our analysis to spatial scales much smaller than the 
horizon size, that is to say when %/k <C 1. Furthermore, 
we consider the weak field approximation, in which the 
scalar fluctuation 5(j> and the gravitational potential $„ 
are considered to be small quantities and only enter the 
equations up to linear order. We consider perturbation 
terms up to second order in v„ (r, x) , as in the Newto- 
nian limit velocities are small with respect to light speed 
(i.e. <C 1 according to our convention). The stress energy 
tensor of the coupled neutrino fluid, T.A, then takes on 
the following form 



= ~Pu (1 + *„) - B{p v ,p v ,5 v W (44) 
T {v) i Q = -B{p v ,p v ,& v ){v v ) i 

r M ° = (i-4$)B(^,^,^)K) 2 

T {y)i = (Pv + S Pv) S i +B(Pv,Pv,6 v )(v 1 ,)i (v„) J 

where we have defined 

B(p u ,p u ,5 v ) = (pu +p„ + p„5 v (l + c 2 )) (45) 

and in the linear regime c 2 = Sp^/Sp^ corresponds to the 
squared sound velocity of neutrino perturbations. 

In the Appendix, perturbation equations are de- 
rived from the conservation equation (TIT)) . V 7 TJ = 
~jjT7 t d a (j). When pressure terms are neglected and for 
Sp c <C 8p v , these equations reduce to the Navier Stokes 
set of equations (fT3]-[T6]), in which the dragging term 
/3<5f>'v„ explicitly appears. The latter was found to be a 
small effect for the chosen values of the coupling, since 
within the considered redshift range the fast oscillations 
of the scalar field average out and, as discussed before, <j> 
can roughly be considered to be constant. 
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VIII. CONCLUSIONS 

In growing neutrino cosmologies, the neutrino mass 
grows in time as a function of a light dark energy scalar 
field, the cosmon. The resulting coupling between cos- 
mon and neutrinos modifies the evolution of the cosmon 
as soon as neutrinos become nonrelativistic. From this 
time on the cosmon evolves only very slowly, and its po- 
tential energy acts similar to a cosmological constant. 
This can provide a natural explanation of the coincidence 
problem without effectively introducing new parameters. 
The current dark energy density can be related to the 
neutrino mass. In addition, the interaction gives rise to 
a new long-range attractive force for neutrinos. This fifth 
force is responsible for a rapid growth of neutrino pertur- 
bations. These perturbations become nonlinear on large 
length scales and may eventually form very large stable 
structures, neutrino lumps. 

We have extended previous studies on the topic to the 
nonlinear regime, providing a method to investigate the 
formation of neutrino lumps by evaluating Navier Stokes 
equations, in which the dragging term due to the cou- 
pling has been suitably included. As we have verified, 
these are the most general equations to compute the 
spacetime evolution of perturbations subject to exter- 
nal forces within the Newtonian limit. Unlike spherical 
collapse methods, the full nonlinear equations describe 
coherently the growth of structures due to an external 
force stronger than gravity. We note that although our 
calculations were performed within a growing neutrino 
scenario, the results are not limited to this case. In- 
stead, the method can be applied for general coupled 
quintessence models, or whenever a force different from 
gravity is present, driving the growth of nonrelativistic 
matter perturbations. 

We have numerically solved the set of hydrodynamical 
equations on a three-dimensional spatial grid. This has 
allowed us to follow the evolution of a single neutrino 
lump in physical space and to determine its properties 
as it approaches stability. Nonrelativistic neutrinos de- 
couple from the background expansion and collapse into 
virialized stable structures. We have identified the char- 
acteristic gravitational potential of these structures as 
a function of redshift and lump size and we have illus- 
trated the detailed behavior of the lumps from the linear 
regime to virialization. Indeed, after a period of fifth- 
force driven collapse, the build-up of large nonradial ve- 
locities leads to stabilization of the lumps. By solving on 
a three-dimensional grid with a general initial profile, we 
are able to provide a detailed illustration of the evolution 
at different distances from the center of the lump, follow- 
ing the change in kinetic and potential energy as well as 
the increasing contribution of nonradial velocities. We 
have also estimated the density profile of the virialized 
neutrino lumps and the associated profile of the gravita- 



tional potential. Limitations of our results for character- 
istic properties of simple neutrino lumps arise from two 
issues. First, at the present stage the numerical precision 
of our algorithm does not allow us to compute the behav- 
ior after virialization sets in and to follow the evolution 
until the neutrino lump becomes approximately stable. 
Second, our work concentrates on the evolution of single 
lumps while we have not addressed interesting topics on 
cosmological scales such as the dynamics of the lumps 
and their possible merging. 

The cosmological gravitational potential provides an 
important mean to test the model versus observations 
and can be used to estimate the impact of neutrino lumps 
on the CMB angular power spectrum. An extrapolation 
from the gravitational potential for single lumps to the 
neutrino lump induced cosmological potential <& v {k) at a 
given wave number k involves an averaging over the dis- 
tribution of neutrino lumps. The present work can only 
give a very rough estimate of the result of this averag- 
ing. Nevertheless, it is apparent that the values of 
resulting from the nonlinear hydrodynamical equations 
remain many orders of magnitude below the extrapo- 
lation of the linearized equations if k ^> 10 -3 Mpc -1 . 
We find a typical value = 10 _6 (fc/fco) -2 , where 

the large uncertainty is reflected by the uncertainty in 
ko that we estimate in the range 8.6 • 10 -3 Mpc -1 < 
fco < 5.6 • 1CU 2 Mpc -1 . This estimate has to be modified 
for k 4.5 • 10 -3 Mpc -1 , where neglected effects beyond 
the Newtonian approximation become relevant, and for 
k £ 1.1 Mpc -1 , where neutrino free streaming prevents 
the formation of lumps. Also, a formation history in- 
volving merging may reduce $,,(£;), especially for small 
k corresponding to very large lumps. 

The characteristic size of <fr v (k) found in this paper is 
of an order of magnitude where it may leave imprints on 
the CMB spectrum at small angular momenta. In par- 
ticular, the late ISW effect could be enhanced, leading to 
stronger correlations between temperature fluctuations in 
the CMB and observed large scale structures in the grav- 
itational potential. (The present observational situation 
hints to an enhancement of this correlation by a factor of 
about 2 as compared to the ACDM model [37[.) Small 
oscillations of the CMB spectrum for angular momenta 
I ;$ 50 are also conceivable. Without a more reliable es- 
timate of z) it seems premature to judge if a given 
growing neutrino model remains compatible with obser- 
vation. We recall in this context that the present work 
concentrates on a particular class of growing neutrino 
models (constant f3) and assumes a large average present 
neutrino mass m„(io) ~ 2 eV. We expect that a smaller 
m v reduces the cosmological effects of neutrino lumps 
since the total neutrino fraction of the energy density f2„ 
gets reduced and the neutrino induced effects set in at 
even more recent cosmological times. 

It may be possible to directly observe large neutrino 
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lumps through their gravitational potential. A possible 
indication would be an observation of structures at large 
length scales where fluctuations in the ACDM model are 
expected to remain linear such that substantial over- or 
underdensities are very rare. The gravitational potential 
of large neutrino lumps may influence a substantial frac- 
tion of space. If we are sitting not too far from a neu- 
trino lump (or even within a neutrino lump) this may 
induce a certain anisotropy of the observed sky on the 
largest scales (i.e. small differences between northern 
and southern hemisphere or similar effects). Finally, the 
rapid recent growth of the neutrino induced gravitational 
potential could lead to an enhancement of the peculiar 
velocities as a characteristic effect of the nongravitational 
interactions [38| . Present observations of the large scale 
bulk flow seem to suggest values substantially larger than 
expected in the ACDM model 3^-42|. Finally, the cos- 
mon may also have a coupling to dark matter, substan- 
tially smaller than the cosmon-neutrino coupling. In this 
case the fifth force effects could also influence the be- 
havior of dark matter. Detection of such effects would 
challenge standard ACDM models, providing a hint for 
interacting cosmologies such as growing neutrino models. 

Let us end with the remark that there is a chance 
that the time variation of the neutrino mass could even 
be detected. Indeed, our understanding of structure 

I 

8' v = 3{%~ >3<f>') (w v - c 2 s ) 8 V - (1 + w v )(l + R V 8 U ) (Vv„ - 3$') - W(l - 3w v + 8 V (1 - 3c 2 )) 



formation and other cosmological features places strong 
upper bounds on the neutrino masses in early cosmology, 
say for z > 5. One may argue about the precise location 
of this bound, but an average neutrino mass m v (z > 5) 
of 0.5 eV would certainly have left a strong imprint on 
cosmology which has not been observed. In consequence, 
if the direct searches for a neutrino mass or the neu- 
trinoless double beta decay indicate a present neutrino 
mass larger than 0.5 eV, this could be interpreted as a 
strong signal in favor of a growing neutrino mass. 
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APPENDIX 
Nonlinear relativistic perturbation equations 

The evolution equations for a neutrino overdensity can 
be derived from conservation of the energy-momentum 
tensor (1451): 



-V$„ + (1 + 8$) 



V8p u 



+ 4(1 + to„)(l + R V 8„)V$ - 2[3V8(f> (1 - 3w v + 8 V (\ - 3c 2 )) 



+ w 2 u 8 u c 2 s +w' v + {l- 3w v ){\ + u>„)(l + R v 5 v ){n - M) + 3«„(1 + c 2 s )(w v - c 2 s )(H - p</>') 



(46) 



_ v$+ l-3^ + Ml r 3c^ v ^ 



1 + 4$ 



v„-V 



(1 + w u )(l + R U 8 V ) 
(2 + c 2 s )v u V8p u - p v c 2 s v v V8 v 



p v {l + w v ){l + R V 8 V ) 



VS Pv 



Pu{l + w v ){l + R v 8 u ) 



3(1 -«£)*' + (ft 



„ l-3^ + i?,^(l-3 C 2) 



1 + Rv8v 



ciV-v v - 



5 v c 2 + w' v 



(1 + w u ){l + R U 8 V ) 



1 + R„8 V 

I 



(38(f)' 



(47) 



where we introduced R v = {1 + c 2 ) / {1 + w v ). The quan- 
tity c 2 — 8p v /8p v equals the squared sound speed of neu- 
trino perturbations on a linear level. 



regime. 

The 0-0-component of Einstein's field equations fulfills 
a 2 a 2 

A$= —Sp+—p(v 2 + 2(1 + 5)$) +3ft$' (48) 



One may verify that (|46|) and (|47|) , if decomposed into 
Fourier modes, reduce to (|39t and (|4Tj) in the linear while the perturbed Klein Gordon equation yields 

I 

8cf>" = a 2 Pp u {S u (l - 3c 2 ) + 2(1 - 3w^)$) - 2a 2 $£/ - a 2 8<j)U A4) + AScf - m.8$ + 40'$' . 

I 



(49) 



If pressure and sound speed are negligible with respect to tional to p u — w u p Ul 
the respective densities, we may omit all terms propor- 



c 2 8p v and their derivatives. 
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Subsequently considering the Newtonian limit yields eqs. 
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